clear all
clc
close all
% 读取数据文件并转换为数组
data = readtable('加噪声加速度set3');  % 读取表格数据
data = table2array(data);  % 将表格转换为数组

% 采样频率
Fs = 500;  % 采样频率为 500 Hz
% 设计低通滤波器，截止频率为 120 Hz
cutoff_frequency = 120;  % 截止频率 120 Hz
data = lowpass(data, cutoff_frequency, 2000);
% 小波包分解
wavelet = 'db4';  % 小波基函数
level = 3;        % 分解层数
wpt = wpdec(data, level, wavelet);

% 获取第三层的节点
layer = 3; % 第三层
tn = leaves(wpt, layer);

% 初始化重构信号
signal_reconstructed = zeros(size(data));

% 计算总信号能量
total_energy = sum(data.^2);

% 能量熵阈值
entropy_threshold = 0.06;  % 设置能量熵的阈值

% 创建日志文件，记录剔除信息
log_file = fopen('node_removal_log.txt', 'w');
fprintf(log_file, 'Node Removal Log\n');
fprintf(log_file, '--------------------------\n');

figure;
sgtitle('Node Signals with Color Indicating Removal Status');

% 提取第三层所有节点，基于能量熵剔除不符合的节点
valid_nodes = [];
for i = 1:length(tn)
    % 获取当前节点的重构信号
    node_signal = wprcoef(wpt, tn(i));      
    % 计算节点的能量
    node_energy = sum(node_signal.^2);
    
    % 计算节点的能量熵
    energy_ratio = node_energy / total_energy;
    entropy = -energy_ratio * log2(energy_ratio + eps);  % 防止对数零错误
    
    % 判断节点是否保留
    is_removed = false;
    if (entropy >= entropy_threshold) && (energy_ratio<0.3)
                fprintf(log_file, 'Node %d: Removed due to energy entropy (%.4f > %.4f)\n', tn(i), entropy, entropy_threshold);
        is_removed = true;
    else
        valid_nodes = [valid_nodes, tn(i)];  % 保留符合条件的节点
    end
    
    % 显示当前节点的信号并使用颜色区分是否剔除
    subplot(length(tn), 1, i);
    plot(node_signal);
    title(['Node ', num2str(tn(i)), ' Signal']);
    xlabel('Time (samples)');
    ylabel('Amplitude');
    
    % 使用颜色区分是否剔除
    if is_removed
        % 使用红色表示剔除的节点
        text(floor(length(node_signal) * 0.7), max(node_signal) * 0.8, 'Removed', 'Color', 'red');
    else
        % 使用绿色表示保留的节点
        text(floor(length(node_signal) * 0.7), max(node_signal) * 0.8, 'Retained', 'Color', 'green');
        % 重构信号
        signal_reconstructed = signal_reconstructed + node_signal;
    end
end




% 初始化全局最大幅值
global_max_amplitude = 0;

% 第一步：遍历所有有效节点，找出所有节点中的最大幅值
for i = 1:length(valid_nodes)
    % 获取当前有效节点的信号
    node_signal = wprcoef(wpt, valid_nodes(i));
    
    % 计算节点频谱
    N_node = length(node_signal);
    U_node = abs(fft(node_signal));
    U_node_half = U_node(1:floor(N_node/2));  % 只取频谱的一半
    
    % 更新全局最大幅值
    node_max_amplitude = max(U_node_half);
    if node_max_amplitude > global_max_amplitude
        global_max_amplitude = node_max_amplitude;
    end
end




% 第二步：根据全局最大幅值对每个节点进行剔除
remaining_nodes = [];  % 初始化存储保留下来的节点

for i = 1:length(valid_nodes)
    % 获取当前有效节点的信号
    node_signal = wprcoef(wpt, valid_nodes(i));
    
    % 计算节点频谱
    N_node = length(node_signal);
    f_node = (0:N_node-1) * (Fs / N_node);
    U_node = abs(fft(node_signal));
    U_node_half = U_node(1:floor(N_node/2));  % 只取频谱的一半
    
    % 找到当前节点的最大幅值
    max_amplitude_node = max(U_node_half);
    
    % 如果当前节点的最大幅值大于全局最大幅值的一半，则剔除整个节点
    if max_amplitude_node > global_max_amplitude / 2
        % 记录剔除节点的日志信息
        fprintf(log_file, 'Node %d: Removed due to max amplitude (%.4f > %.4f)\n', ...
                valid_nodes(i), max_amplitude_node, global_max_amplitude / 2);
    else
        % 如果节点未被剔除，则保留节点并加入重构信号
        remaining_nodes = [remaining_nodes, valid_nodes(i)];  % 保留该节点
        signal_reconstructed = signal_reconstructed + node_signal;  % 累加重构信号
    end
end

% 绘制保留节点的频谱
figure;
for i = 1:length(remaining_nodes)
    node_signal = wprcoef(wpt, remaining_nodes(i));
    N_node = length(node_signal);
    f_node = (0:N_node-1) * (Fs / N_node);
    U_node = abs(fft(node_signal));
    U_node_half = U_node(1:floor(N_node/2));  % 只取频谱的一半
    f_half = f_node(1:floor(N_node/2));       % 对应的频率轴
    
    % 绘制保留节点的频谱
    subplot(length(remaining_nodes), 1, i);
    plot(f_half, U_node_half);
    title(['Node ', num2str(remaining_nodes(i)), ' Frequency Spectrum']);
    xlabel('Frequency (Hz)');
    ylabel('Amplitude');
end

% 关闭日志文件
fclose(log_file);
% 绘制保留节点的频谱
figure;
for i = 1:length(remaining_nodes)
    node_signal = wprcoef(wpt, remaining_nodes(i));
    N_node = length(node_signal);
    f_node = (0:N_node-1) * (Fs / N_node);
    U_node = abs(fft(node_signal));
    U_node_half = U_node(1:floor(N_node/2));  % 只取频谱的一半
    f_half = f_node(1:floor(N_node/2));       % 对应的频率轴
    
    % 绘制保留节点的频谱
    subplot(length(remaining_nodes), 1, i);
    plot(f_half, U_node_half);
    title(['Node ', num2str(remaining_nodes(i)), ' Frequency Spectrum']);
    xlabel('Frequency (Hz)');
    ylabel('Amplitude');
end
% 显示原始信号和降噪后的信号
figure;
subplot(2,1,1);
plot(data);
title('Original Signal');
xlabel('Time (samples)');
ylabel('Amplitude');

subplot(2,1,2);
plot(signal_reconstructed);
title('Denoised Signal after Frequency Amplitude Processing');
xlabel('Time (samples)');
ylabel('Amplitude');

% 设计低通滤波器，截止频率为 120 Hz
cutoff_frequency = 120;  % 截止频率 120 Hz
signal_reconstructed = lowpass(signal_reconstructed, cutoff_frequency, Fs);



% 设置相同的Y轴范围
max_value = max([max(data), max(signal_reconstructed)]);
min_value = min([min(data), min(signal_reconstructed)]);
if min_value < max_value
    subplot(2,1,1);
    ylim([min_value, max_value]);

    subplot(2,1,2);
    ylim([min_value, max_value]);
else
    warning('min_value should be less than max_value.');
end

% 变分模态分解（VMD）

% 参数设置
alpha = 1760;       % 路由参数，调节每个模态的稀疏性
tau = 0;            % 噪声容限
K = 30;              % 模态数
DC = 0;             % 不对第一个模态进行DC偏置
init = 1;           % 初始化方式，1表示从均匀分布开始
tol = 1e-7;         % 收敛容限

% 使用VMD函数进行分解
[u, u_hat, omega] = VMD(signal_reconstructed, alpha, tau, K, DC, init, tol);

% 显示分解后的模态以及原始信号
global_min = min([u(:); min(signal_reconstructed)]);  % 包含原始信号的最小值
global_max = max([u(:); max(signal_reconstructed)]);  % 包含原始信号的最大值

figure;

% 绘制原始信号
subplot(K+1, 1, 1);
plot(signal_reconstructed);
title('Reconstructed Signal');
xlabel('Time (samples)');
ylabel('Amplitude');
ylim([global_min, global_max]);  % 统一Y轴范围

% 绘制分解后的模态
for i = 1:K
    subplot(K+1, 1, i+1);
    plot(u(i, :));
    title(['IMF', num2str(i)]);
    xlabel('Time (samples)');
    ylabel(['IMF', num2str(i)]);
    ylim([global_min, global_max]);  % 统一Y轴范围
end

% 计算并显示各模态的频域图
N = length(signal_reconstructed);  % 获取信号长度
f = (0:N-1)*(Fs/N);  % 计算频率轴
figure;

% 初始化存储每个模态中心频率的数组
center_frequencies_peak = zeros(K, 1);

% 绘制原始信号的频谱
subplot(K+1, 1, 1);
U_original = abs(fft(signal_reconstructed));
U_original = U_original(1:floor(N/2));
f_half = f(1:floor(N/2));
plot(f_half, U_original);
title('Frequency Spectrum of Reconstructed Signal');
xlabel('Frequency (Hz)');
ylabel('Amplitude');

% 绘制每个模态的频谱并找出中心频率（最高峰值）
for i = 1:K
    % 计算傅里叶变换
    U = abs(fft(u(i, :)));
    
    % 只显示一半频率范围（0 到 Fs/2）
    U_half = U(1:floor(N/2));  % 取频率的前半部分
    f_half = f(1:floor(N/2));  % 对应的频率轴
    
    % 找到频谱中的最高峰值及其对应的频率
    [peak_value, peak_idx] = max(U_half);  % 获取峰值和峰值位置
    center_frequencies_peak(i) = f_half(peak_idx);  % 记录最高峰对应的频率 
end

% 按中心频率升序排序IMF
[center_frequencies_peak_sorted, idx] = sort(center_frequencies_peak);
u_sorted = u(idx, :);  % 对IMF按中心频率排序

% 绘制每个模态的频谱并找出中心频率（最高峰值）
for i = 1:K
    % 绘制排序后的IMF频谱
    U = abs(fft(u_sorted(i, :)));
    U_half = U(1:floor(N/2));
    f_half = f(1:floor(N/2));
    
    % 绘制频域图
    subplot(K+1, 1, i+1);
    plot(f_half, U_half);
    title(['Frequency Spectrum of IMF ', num2str(i)]);
    xlabel('Frequency (Hz)');
    ylabel(['Amplitude of IMF ', num2str(i)]);
    
    % 标记出峰值位置
    hold on;
    [peak_value, peak_idx] = max(U_half);
    plot(f_half(peak_idx), peak_value, 'ro');  % 用红点标出峰值
    text(f_half(peak_idx), peak_value, sprintf('Peak: %.2f Hz', f_half(peak_idx)), ...
        'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');
end

% 输出每个排序后的IMF的中心频率
fprintf('Sorted Center frequencies (peak values):\n');
for i = 1:K
    fprintf('Sorted IMF %d: Peak Frequency = %.2f Hz\n', i, center_frequencies_peak_sorted(i));
end

% 绘制IMF序号与中心频率的关系图
figure;
plot(1:K, center_frequencies_peak_sorted, '-o', 'LineWidth', 2, 'MarkerSize', 8);
title('Center Frequency vs. IMF Number (Sorted)');
xlabel('IMF Number');
ylabel('Center Frequency (Hz)');
grid on;

